En este trabajo, se describe una aplicación práctica de las Máquinas de Soporte Vectorial, como soporte para acompañar el desarrollo teórico realizado en el primer video de la serie. Se exploran alternativas de funciones kernel, impacto de los parámetros y rendimiento de clasificadores en el caso de estudio: identificación prematura de pacientes con enfermedad PDAC. Se comparan los algoritmos con un modelo de regresión logstica.
El cáncer pancreático es una neoplasia extremadamente mortal, con una incidencia del 30% de los tumores diagnosticados a nivel global durante el año 2020. Una vez diagnosticado, la tasa de supervivencia en los siguientes 5 años es de apenas el 10%. No obstante, si el mismo es detectado en etapas tempranas, la probabilidad de supervivencia aumenta. Desafortunadamente, muchos casos de cáncer pancreático no muestran síntomas detectables sino hasta que el mismo ha generado metástasis. En este sentido, un diagnóstico temprano de los pacientes con cáncer pancreático representa una necesidad clínica insatisfecha y una ventana de oportunidad para el desarrollo de nuevas herramientas que asistan a la identificación temprana de esta afección.
Por tal motivo, se eligió aplicar una clasificación supervisada con Máquinas de Soporte vectorial sobre un conjunto de datos publicado recientemente (Debernardi y col, 2020), con el objeto de ilustrar la fortaleza de dicho algoritmo para clasificar pacientes en función de diversos marcadores urinarios, previo al diagnóstico de PDAC.
La base de datos de Debernardi y colaboradores (doi: 10.1371/journal.pmed.1003489.) cuenta con la edad (años), sexo, valores de los biomarcadores urinarios creatinina [mg/ml], LYVE1 [ng/ml], REG1B [ng/ml], y TFF1 [ng/ml] para pacientes en tres categorías de la variable diagnóstico: sanos, con alguna condición benigna, y con PDAC (adenocarcinoma pancreático ductal, el tipo de cáncer pancreático más recurrente). Los datos fueron generados en 4 centros de investigación. En este trabajo se decidió utilizar los datos de pacientes cuyo diagnóstico era “normal” (sanos) o “maligno” (PDAC), conformando una base de datos de 381 registros, y se utilizó la totalidad de los datos sin tener en cuenta su proveniencia.
Fuente original de los datos: https://doi.org/10.1371/journal.pmed.1003489.s009
knitr::opts_chunk$set(echo = TRUE)
library(mlr)
## Loading required package: ParamHelpers
## Warning message: 'mlr' is in 'maintenance-only' mode since July 2019.
## Future development will only happen in 'mlr3'
## (<https://mlr3.mlr-org.com>). Due to the focus on 'mlr3' there might be
## uncaught bugs meanwhile in {mlr} - please consider switching.
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
#remotes::install_github("vqv/ggbiplot")
library(ggbiplot)
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## Loading required package: scales
## Loading required package: grid
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(ggpubr)
##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
##
## mutate
library(knitr)
library(googlesheets)
library(tidyr)
library(plotly)
##
## Attaching package: 'plotly'
## The following objects are masked from 'package:plyr':
##
## arrange, mutate, rename, summarise
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
theme <- theme(text = element_text(size=10),plot.title = element_text(size=12, face="bold.italic",
hjust = 0.5), axis.title.x = element_text(size=10, face="bold", colour='black'),
axis.title.y = element_text(size=10, face="bold"),panel.border = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.title = element_text(face="bold"))
id <- '11wXwKESJrMw6_hqljbOzpn7ncXPjNWeV' # nombre del archivo en google drive
data <- read.csv(sprintf("https://docs.google.com/uc?id=%s&export=download", id))
#data <- read.csv('pancreas.csv')
# Cambio el nombre de las variables a español
colnames(data) <- c('id','cohorte','origen','edad','sexo','diagnosis','estadio','diagnosis_benigno','plasma_CA19_9','creatinina','LYVE1','REG1B','TFF1','REG1A')
# Elijo normal y maligno, y renombro 1=normal, 3=maligno
data <- data %>% filter ( diagnosis=='1' | diagnosis=='3')
data <- data %>% mutate(diagnosis = sub("1", "normal", diagnosis))
data <- data %>% mutate(diagnosis = sub("3", "maligno", diagnosis))
# Transformo como factor
data$sexo <- as.factor(data$sexo)
data$diagnosis <- factor(data$diagnosis, levels= c('normal','maligno'))
# Me quedo sólo con columnas: edad,sexo,diagnosis,creatinina,LYVE1,REG1B,TFF1 y las reordeno
data <- data %>% dplyr::select(6,5,4,10:13)
# Normalizo variables
#data <- data %>% mutate_at(c('edad', 'creatinina','LYVE1','REG1B','TFF1'), ~(scale(.) %>% as.vector))
# Imprimo tabla
kable(data.frame(variable = names(data),
class = sapply(data, class),
primeros_valores = sapply(data, function(x) paste0(head(x), collapse = ", ")),
row.names = NULL))
| variable | class | primeros_valores |
|---|---|---|
| diagnosis | factor | normal, normal, normal, normal, normal, normal |
| sexo | factor | F, F, M, M, M, M |
| edad | integer | 33, 81, 51, 61, 62, 53 |
| creatinina | numeric | 1.83222, 0.97266, 0.78039, 0.70122, 0.21489, 0.84825 |
| LYVE1 | numeric | 0.8932192, 2.037585, 0.1455889, 0.00280488, 0.00085956, 0.003393 |
| REG1B | numeric | 52.94884, 94.46703, 102.366, 60.579, 65.54, 62.126 |
| TFF1 | numeric | 654.282174, 209.48825, 461.141, 142.95, 41.088, 59.793 |
# Gráfico ggpairs
data%>% ggpairs(.,mapping=ggplot2::aes(color = diagnosis,alpha = 0.1),
upper = list(continuous = wrap("cor", size = 2.5),discrete = "blank", combo="blank"),
lower = list(combo = "box"),progress = F)+
theme+
labs(title= 'Descripción de variables en la base de datos', x='Variable', y='Variable')+
scale_fill_manual(values=c('royalblue2','red'))+
scale_color_manual(values=c('royalblue2','#ff7474ff'))
set.seed(1409) # para asegurar reproducibilidad
dt = sort(sample(nrow(data), nrow(data)*.7))
datos_tr<-data[dt,]
datos_te<-data[-dt,]
Esta función toma como argumentos el modelo, los datos de entrenamiento y prueba. Calcula y muestra métricas, y curvas ROC
metricas <- function(modelo, train, test, nombre){
pred_svm_test= predict(modelo, newdata = test)
acc_svm_test <- round(measureACC(as.data.frame(pred_svm_test)$truth, as.data.frame(pred_svm_test)$response),3)
AUC_svm_test <- round(measureAUC(as.data.frame(pred_svm_test)$prob.maligno,as.data.frame(pred_svm_test)$truth,'normal','maligno'),3)
# ···························································································
# Predicción TRAIN (naive)
pred_svm_train = predict(modelo, newdata = train) # por si quiero ver naive sobre training
acc_svm_train <- round(measureACC(as.data.frame(pred_svm_train)$truth, as.data.frame(pred_svm_train)$response),3)
AUC_svm_train <- round(measureAUC(as.data.frame(pred_svm_train)$prob.maligno, as.data.frame(pred_svm_train)$truth, 'normal','maligno'),3)
# ···························································································
# Cambio el threshold [esto lo hago para train y test]
acc=NULL
acc2=NULL
threshold = seq(0.1,0.95,0.01)
for (i in 1:length(threshold)) {
pred = setThreshold(pred_svm_test, threshold = threshold[i])
acc[i] = measureACC(as.data.frame(pred)$truth, as.data.frame(pred)$response)}
for (i in 1:length(threshold)) {
pred2 = setThreshold(pred_svm_train, threshold = threshold[i])
acc2[i] = measureACC(as.data.frame(pred2)$truth, as.data.frame(pred2)$response)}
par(mfcol = c(1,2))
new_df1 <- as.data.frame(cbind(threshold,acc))
new_df1 <- new_df1%>%mutate(sub_data='test')
new_df2 <- as.data.frame(cbind(threshold,acc2))
colnames(new_df2) <- c('threshold','acc')
new_df2 <- new_df2%>%mutate(sub_data='train')
new_df <- as.data.frame(rbind(new_df1,new_df2))
#new_df1[which.max(new_df1$acc),"threshold"] # test 0.39
#new_df2[which.max(new_df2$acc),"threshold"] # train 0.53
# ···························································································
# Gráfico de cómo varía la métrica de performance accuracy, de acuerdo al umbral elegido
plot_acc <- ggplot(new_df, aes(x=threshold, y=acc)) + geom_line(aes(color = sub_data,linetype=sub_data)) +
theme +labs(x='Umbral', y='Métrica de performance (accuracy)',
title= 'Evaluación del modelo de Máquinas de soporte vectorial SVM') +
scale_color_manual(values = c("red", "darkred"),labels=c('prueba','entrenamiento')) +
scale_linetype_manual(values=c(1,2), labels=c('prueba','entrenamiento')) +
labs(color='Conjunto de\n evaluación',linetype='Conjunto de\n evaluación')
print(plot_acc)
# Para independizarnos de la elección del umbral, grafico curvas ROC para las predicciones del modelo SVM con los datos de TEST y TRAIN
df_svm = generateThreshVsPerfData(list(svm_te = pred_svm_test, svm_tr = pred_svm_train),
measures = list(fpr, tpr, mmce))
plot_roc <- plotROCCurves(df_svm) + theme +
labs(title=paste0('Curva ROC del modelo - ', nombre),
x='Tasa de falsos positivos (FPR)', y='Tasa de positivos verdaderos (TPR)',
color='Conjunto de\n evaluación') +
scale_color_manual(values = c("red", "darkred"), labels=c('prueba','entrenamiento'))
# geom_label(label="AUC= 0.894", x=0.35, y=0.75, label.size = 0.3, size=4,
# color = "red",fill="white") +
# geom_label(label="AUC= 0.925", x=0.07, y=0.97, label.size = 0.3, size=4,
# color = "darkred",fill="white")
print(plot_roc)
# ················ Métricas del modelo de SVM ················
Métrica <- c('valor','datos')
Accuracy <- c(acc_svm_test,'prueba')
Accuracy. <- c(acc_svm_train,'entrenamiento')
AUC_ROC <- c(AUC_svm_test,'prueba')
AUC_ROC. <- c(AUC_svm_train,'entrenamiento')
# Imprimo resultados
kable(rbind(Métrica, Accuracy, Accuracy., AUC_ROC, AUC_ROC.))
}
Se utilizará la librería mlr, que ofrece herramientas para trabajar con experimentos de machine learning. En el caso de svm, mlr utiliza las funciones del paquete e1071.
Daremos un primer vistazo a un modelo simple con kernel lineal. Para ello se definen la tarea de clasificación y el modelo a ajustar. Posteriormente, se entrena el modelo con los datos previamente reservados par ael entrenamiento.
# Defino modelo SVM
set.seed(1)
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis")
lrn_svmL = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "linear", cost=1))
mod_svmL = mlr::train(lrn_svmL, task)
Calcularemos nuestras primeras métricas sobre este modelo, para conocer el desempeño
metricas(mod_svmL, datos_tr[-2], datos_te[-2], 'SVM con Kernel Lineal')
| Métrica | valor | datos |
| Accuracy | 0.826 | prueba |
| Accuracy. | 0.872 | entrenamiento |
| AUC_ROC | 0.891 | prueba |
| AUC_ROC. | 0.924 | entrenamiento |
¿Cómo se ve la frontera de decisión para este clasificador? Dada la dimensionalidad de nuestro dataset, no es posible observarla de manera directa. A modo ilustrativo, haremos una reducción de la dimensionalidad a través del cómputo de las Componentes Principales.
# Preparo los datos para análisis de componentes principales
datos_para_acp = data[c(3:7)] # todas las variables numéricas
datos.pc = prcomp(datos_para_acp,scale = TRUE) #escalo los datos
summary(datos.pc)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.6043 1.0609 0.7881 0.63208 0.52941
## Proportion of Variance 0.5147 0.2251 0.1242 0.07991 0.05605
## Cumulative Proportion 0.5147 0.7398 0.8640 0.94395 1.00000
Conservaremos sólo las dos primeras componentes, con relativa tranquilidad de capturar con ellas casi el 74% de la varianza de los datos -más que suficiente para un análisis ilustrativo.
data_toplot = datos.pc$x[,c('PC1','PC2')]
data_toplot = data.frame(data_toplot, data[,'diagnosis'])
colnames(data_toplot) = c('PC1','PC2', 'diagnosis')
Definiremos una pequeña función para crear una grilla de datos. Con ella, generaremos datos espaciados uniformemente y cubriendo la totalidad del rango de nuestras variables en la muestra estudiada. Luego, alimentaremos nuestro clasificador con la grilla, y los puntos servirán de testigo para dibujar una frontera de decisión.
make.grid = function(x, n = 75) {
grange = apply(x, 2, range)
x1 = seq(from = grange[1,1], to = grange[2,1], length = n)
x2 = seq(from = grange[1,2], to = grange[2,2], length = n)
expand.grid(X1 = x1, X2 = x2)
}
set.seed(1)
task = makeClassifTask(data = data_toplot, target = "diagnosis")
lrn_plot = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "linear", cost = 1))
mod_plot = mlr::train(lrn_plot, task)
xgrid = make.grid(data_toplot[, c('PC1','PC2')])
ygrid = predict(mod_plot, newdata = xgrid)
plot(xgrid, col = c("red","green")[as.numeric(unlist(ygrid))], pch = 20, cex = .2)
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
points(data_toplot[,c('PC1','PC2')], col = as.integer(data_toplot$diagnosis) + 1 , pch = 19)
# Con esta línea, se pueden conocer los vectores de soporte
#points(data_toplot[,c('PC1','PC2')][mod_plot$learner.model$index,], pch = 5, cex = 2)
En esta figura, se puede ver el diagnóstico real para cada paciente. Y, además, gracias a nuestra grilla podemos ver la clase que el algoritmo asignaría a cualquier otro hipotético punto que quisiéramos incorporar. Los datos están relativamente condensados pero muy superpuestos. Evidentemente, un análisis bidimensional no es suficiente para separarlos -al menos no el obtenido con esta metodología.
Los kernels que ofrece e1071, se pueden conocer a través de la función getParamSet. Además también, si quisiéramos, podríamos investigar los parámetros disponibles para acompañar a cada kernel particular, a través de este mismo objeto.
getParamSet("classif.svm")$pars$kernel$values
## $linear
## [1] "linear"
##
## $polynomial
## [1] "polynomial"
##
## $radial
## [1] "radial"
##
## $sigmoid
## [1] "sigmoid"
¿Cómo se ven las fronteras de separación con otros kernels más flexibles? Si bien no podemos esperar una separación perfecta en dos dimensiones, quizás haya vectores cuyas clases reales sean mejor capturadas por otros kernels.
En primer lugar, observaremos kernels polinomiales de grados entre 2 y 5. Ya puede verse una flexibilidad tremendamente superior al elemental kernel lineal aunque hay algunos puntos difíciles de capturar.
set.seed(1)
for(i in 2:5){
task = makeClassifTask(data = data_toplot, target = "diagnosis")
lrn_plot = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "polynomial", degree=i, cost = 1))
mod_plot = mlr::train(lrn_plot, task)
xgrid = make.grid(data_toplot[, c('PC1','PC2')])
ygrid = predict(mod_plot, newdata = xgrid)
plot(xgrid, col = c("red","green")[as.numeric(unlist(ygrid))], pch = 20, cex = .2)
points(data_toplot[,c('PC1','PC2')], col = as.integer(data_toplot$diagnosis) + 1 , pch = 19)
}
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
¿Cuál es el efecto del parámetro C? Observaremos un barrido entre unos pocos valores del parámetro, para un kernel polinomial de grado 4. Marca la diferencia entre un ajuste más suave o más estricto.
set.seed(1)
for(i in 2:5){
task = makeClassifTask(data = data_toplot, target = "diagnosis")
lrn_plot = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "polynomial", degree=4, cost = i))
mod_plot = mlr::train(lrn_plot, task)
xgrid = make.grid(data_toplot[, c('PC1','PC2')])
ygrid = predict(mod_plot, newdata = xgrid)
plot(xgrid, col = c("red","green")[as.numeric(unlist(ygrid))], pch = 20, cex = .2)
points(data_toplot[,c('PC1','PC2')], col = as.integer(data_toplot$diagnosis) + 1 , pch = 19)
}
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
A continuación, observaremos las fronteras que describe el kernel Radial. Aprovecharemos, simultáneamente, para tener una noción gráfica del efecto del parámetro gamma. Si recordamos la ecuación que define la función del kernel radial, podemos esperar una influencia menor de la distancia entre dos puntos cuanto más grande es gamma
Puede observarse todavía más flexibilidad que en el caso polinomial, y se ve claramente cómo varia la frontera ante el barrido de gamma. Se muestran valores de gamma entre 1 y 6, para un valor constante de C. El efecto de localidad es mayor cada vez.
for(i in 1:6){
set.seed(1)
task = makeClassifTask(data = data_toplot, target = "diagnosis")
lrn_plot = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "radial", cost=0.8, gamma = i))
mod_plot = mlr::train(lrn_plot, task)
xgrid = make.grid(data_toplot[, c('PC1','PC2')])
ygrid = predict(mod_plot, newdata = xgrid)
plot(xgrid, col = c("red","green")[as.numeric(unlist(ygrid))], pch = 20, cex = .2)
points(data_toplot[,c('PC1','PC2')], col = as.integer(data_toplot$diagnosis) + 1 , pch = 19)
}
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
#points(datos_tr[,c(x1,x2)][mod_plot$learner.model$index,], pch = 5, cex = 2)
Finalmente, observaremos la frontera del kernel Sigmoide. La curva que dibuja es altamente flexible, pero no ha capturado alunos puntos específicos ¿Será posible obtener un mejor ajuste?
set.seed(1)
task = makeClassifTask(data = data_toplot, target = "diagnosis")
lrn_plot = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "sigmoid", cost=2, gamma = 0.5))
mod_plot = mlr::train(lrn_plot, task)
xgrid = make.grid(data_toplot[, c('PC1','PC2')])
ygrid = predict(mod_plot, newdata = xgrid)
plot(xgrid, col = c("red","green")[as.numeric(unlist(ygrid))], pch = 20, cex = .2)
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
points(data_toplot[,c('PC1','PC2')], col = as.integer(data_toplot$diagnosis) + 1 , pch = 19)
#points(datos_tr[,c(x1,x2)][mod_plot$learner.model$index,], pch = 5, cex = 2)
La definición de C, gamma, el grado -para polinomiales- e incluso el kernel óptimo a utilizar puede hacerse a través de una búsqueda por grilla, aleatoria o bayesiana. A modo de ejemplo, buscaremos el set óptimo de parámetros con una búsqueda aleatoria y cross validation.
El paquete mlr nos sirve en este caso para ejecutar el flujo de trabajo de forma compacta y ordenada. Definida la tarea y el clasificador, definimos el espacio de búsqueda paramétrica. Luego, definimos la estrategia de búsqueda y cv.
# Tarea y clasificador
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis")
svm <- makeLearner("classif.svm")
# Parámetros a optimizar
kernels <- c("polynomial", "radial", "sigmoid")
svmParamSpace <- makeParamSet(
makeDiscreteParam("kernel", values = kernels),
makeIntegerParam("degree", lower = 1, upper = 3),
makeNumericParam("cost", lower = 0.1, upper = 10),
makeNumericParam("gamma", lower = 0.1, 10))
# Esquema de búsqueda
randSearch <- makeTuneControlRandom(maxit = 30)
cvForTuning <- makeResampleDesc("Holdout", split = 2/3)
#Parámetros óptimos
tunedSvmPars <- tuneParams("classif.svm", task = task,
resampling = cvForTuning,
par.set = svmParamSpace,
control = randSearch)
## [Tune] Started tuning learner classif.svm for parameter set:
## Type len Def Constr Req Tunable Trafo
## kernel discrete - - polynomial,radial,sigmoid - TRUE -
## degree integer - - 1 to 3 - TRUE -
## cost numeric - - 0.1 to 10 - TRUE -
## gamma numeric - - 0.1 to 10 - TRUE -
## With control class: TuneControlRandom
## Imputation value: 1
## [Tune-x] 1: kernel=sigmoid; degree=3; cost=1.27; gamma=0.487
## [Tune-y] 1: mmce.test.mean=0.2471910; time: 0.0 min
## [Tune-x] 2: kernel=polynomial; degree=2; cost=9.45; gamma=5.17
## [Tune-y] 2: mmce.test.mean=0.2359551; time: 0.0 min
## [Tune-x] 3: kernel=radial; degree=3; cost=4.07; gamma=8.97
## [Tune-y] 3: mmce.test.mean=0.1910112; time: 0.0 min
## [Tune-x] 4: kernel=radial; degree=1; cost=1.41; gamma=2.73
## [Tune-y] 4: mmce.test.mean=0.1235955; time: 0.0 min
## [Tune-x] 5: kernel=radial; degree=1; cost=8.45; gamma=1.75
## [Tune-y] 5: mmce.test.mean=0.1573034; time: 0.0 min
## [Tune-x] 6: kernel=sigmoid; degree=2; cost=2.8; gamma=3.5
## [Tune-y] 6: mmce.test.mean=0.3483146; time: 0.0 min
## [Tune-x] 7: kernel=sigmoid; degree=2; cost=3.43; gamma=9.68
## [Tune-y] 7: mmce.test.mean=0.3707865; time: 0.0 min
## [Tune-x] 8: kernel=radial; degree=1; cost=3.64; gamma=5.18
## [Tune-y] 8: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 9: kernel=polynomial; degree=1; cost=2.38; gamma=1.71
## [Tune-y] 9: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 10: kernel=radial; degree=2; cost=1.96; gamma=4.14
## [Tune-y] 10: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 11: kernel=radial; degree=2; cost=9.4; gamma=0.556
## [Tune-y] 11: mmce.test.mean=0.1011236; time: 0.0 min
## [Tune-x] 12: kernel=polynomial; degree=2; cost=9.07; gamma=8.82
## [Tune-y] 12: mmce.test.mean=0.2247191; time: 0.0 min
## [Tune-x] 13: kernel=sigmoid; degree=3; cost=4.09; gamma=5.61
## [Tune-y] 13: mmce.test.mean=0.3370787; time: 0.0 min
## [Tune-x] 14: kernel=radial; degree=1; cost=9.26; gamma=8.8
## [Tune-y] 14: mmce.test.mean=0.1910112; time: 0.0 min
## [Tune-x] 15: kernel=polynomial; degree=3; cost=3.32; gamma=1.69
## [Tune-y] 15: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 16: kernel=polynomial; degree=2; cost=2.29; gamma=0.12
## [Tune-y] 16: mmce.test.mean=0.3820225; time: 0.0 min
## [Tune-x] 17: kernel=sigmoid; degree=3; cost=5.11; gamma=7.29
## [Tune-y] 17: mmce.test.mean=0.3595506; time: 0.0 min
## [Tune-x] 18: kernel=radial; degree=2; cost=3.93; gamma=3.44
## [Tune-y] 18: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 19: kernel=radial; degree=1; cost=6.99; gamma=8.17
## [Tune-y] 19: mmce.test.mean=0.1910112; time: 0.0 min
## [Tune-x] 20: kernel=sigmoid; degree=1; cost=9.07; gamma=4.19
## [Tune-y] 20: mmce.test.mean=0.3483146; time: 0.0 min
## [Tune-x] 21: kernel=polynomial; degree=3; cost=2.49; gamma=3.07
## [Tune-y] 21: mmce.test.mean=0.1685393; time: 0.0 min
## [Tune-x] 22: kernel=polynomial; degree=3; cost=0.237; gamma=5.56
## [Tune-y] 22: mmce.test.mean=0.1685393; time: 0.0 min
## [Tune-x] 23: kernel=radial; degree=3; cost=2.61; gamma=4.17
## [Tune-y] 23: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 24: kernel=radial; degree=1; cost=6.44; gamma=3.53
## [Tune-y] 24: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 25: kernel=radial; degree=3; cost=3.9; gamma=2.36
## [Tune-y] 25: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 26: kernel=sigmoid; degree=1; cost=6.89; gamma=7.94
## [Tune-y] 26: mmce.test.mean=0.3707865; time: 0.0 min
## [Tune-x] 27: kernel=sigmoid; degree=3; cost=3.12; gamma=7.77
## [Tune-y] 27: mmce.test.mean=0.3707865; time: 0.0 min
## [Tune-x] 28: kernel=polynomial; degree=2; cost=8.35; gamma=7.61
## [Tune-y] 28: mmce.test.mean=0.2247191; time: 0.0 min
## [Tune-x] 29: kernel=polynomial; degree=1; cost=3.37; gamma=8.65
## [Tune-y] 29: mmce.test.mean=0.1573034; time: 0.0 min
## [Tune-x] 30: kernel=sigmoid; degree=2; cost=4.42; gamma=1.85
## [Tune-y] 30: mmce.test.mean=0.3370787; time: 0.0 min
## [Tune] Result: kernel=radial; degree=2; cost=9.4; gamma=0.556 : mmce.test.mean=0.1011236
La salida nos indica la mejor combinación, para esta estrategia y espacio de búsqueda con 30 iteraciones.
[Tune] Result: kernel=radial; degree=1; cost=3.37; gamma=2.96 : mmce.test.mean=0.0786517
A continuación pondremos a prueba el desempeño de varios clasificadores, incluído el que reúne los parámetros óptimos. Contrastaremos su desempeño con el de una regresión logística. Ahora, sobre el problema original y no su versión reducida.
Calcularemos, para cada modelo: 1. Curvas de Accuracy para diferentes thresholds, en sets de entrenamiento y prueba. 2. Curvas ROC en los mismos sets. AUC será nuestra medida de mayor interés.
# Defino modelo SVM
set.seed(1)
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis")
lrn_svmL = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "linear", cost=1))
mod_svmL = mlr::train(lrn_svmL, task)
metricas(mod_svmL, datos_tr[-2], datos_te[-2], 'SVM con Kernel Lineal')
| Métrica | valor | datos |
| Accuracy | 0.826 | prueba |
| Accuracy. | 0.872 | entrenamiento |
| AUC_ROC | 0.891 | prueba |
| AUC_ROC. | 0.924 | entrenamiento |
Kernel Lineal: AUC_ROC 0.891 prueba
set.seed(1)
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis")
lrn_svmP = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "polynomial", cost=2 , degree=3))
mod_svmP = mlr::train(lrn_svmP, task)
metricas(mod_svmP, datos_tr[-2], datos_te[-2], 'SVM con Kernel Polinomial (grado 3)')
| Métrica | valor | datos |
| Accuracy | 0.783 | prueba |
| Accuracy. | 0.861 | entrenamiento |
| AUC_ROC | 0.895 | prueba |
| AUC_ROC. | 0.933 | entrenamiento |
Kernel Polinomial: AUC_ROC 0.895 prueba
Parámetros óptimos: [Tune] Result: kernel=radial; degree=1; cost=3.37; gamma=2.96
# Defino modelo SVM
set.seed(1)
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis")
lrn_svmR = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "radial", cost=3.37, gamma =2.96))
mod_svmR = mlr::train(lrn_svmR, task)
metricas(mod_svmR, datos_tr[-2], datos_te[-2], 'Radial')
| Métrica | valor | datos |
| Accuracy | 0.809 | prueba |
| Accuracy. | 0.985 | entrenamiento |
| AUC_ROC | 0.902 | prueba |
| AUC_ROC. | 1 | entrenamiento |
# chequeo el balance de las distintas clases de la variable diagnóstico en el conjunto de todos los datos, los datos de prueba y los datos de entramiento.
Entrenamiento <- table(datos_tr$diagnosis) #126 normal, 140 maligno
Prueba <- table(datos_te$diagnosis) # 57 normal, 58 maligno
Total <- table(data$diagnosis) # 183 normal, 198 maligno
kable(rbind(Entrenamiento, Prueba,Total))
| normal | maligno | |
|---|---|---|
| Entrenamiento | 126 | 140 |
| Prueba | 57 | 58 |
| Total | 183 | 198 |
# Armo modelo de regresión logistica.
set.seed(1)
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis")
lrn = makeLearner("classif.logreg", predict.type = "prob")
mod_lr = mlr::train(lrn, task)
metricas(mod_lr, datos_tr[-2], datos_te[-2], 'Logística')
| Métrica | valor | datos |
| Accuracy | 0.817 | prueba |
| Accuracy. | 0.872 | entrenamiento |
| AUC_ROC | 0.889 | prueba |
| AUC_ROC. | 0.925 | entrenamiento |
kable(mod_lr$learner.model$coefficients)
| x | |
|---|---|
| (Intercept) | -5.4691572 |
| edad | 0.0692457 |
| creatinina | -1.0962220 |
| LYVE1 | 0.4330571 |
| REG1B | 0.0027155 |
| TFF1 | 0.0020305 |
Finalmente, compararemos las métricas calculadas sobre el desempeño de todos los modelos construídos.
# ······················ curvas ROC para todos los modelos ········································
pred_lr = predict(mod_lr, newdata = datos_te[-2])
pred_svmL = predict(mod_svmL, newdata = datos_te[-2])
pred_svmP = predict(mod_svmP, newdata = datos_te[-2])
pred_svmR = predict(mod_svmR, newdata = datos_te[-2])
df_todos = generateThreshVsPerfData(list(lg=pred_lr, svmL = pred_svmL, svmP = pred_svmP, svmR = pred_svmR), measures = list(fpr, tpr, mmce))
plotROCCurves(df_todos) + theme +
labs(title='Curvas ROC de modelos de clasificación supervisada (datos de prueba)',
x='Tasa de falsos positivos (FPR)', y='Tasa de positivos verdaderos (TPR)',
color=' Modelo en\n evaluación') +
scale_color_manual(values = c("red", "black", "blue", "darkgreen"),
labels=c('Reg log','SVM (lineal)','SVM (polinomial)','SVM (radial)'))+
theme(legend.position=c(0.915,0.25))
# ············ valores AUC para todos los modelos cuando se consideran todas las variables ···················
AUC_lg_te <- round(measureAUC(as.data.frame(pred_lr)$prob.maligno, as.data.frame(pred_lr)$truth, 'normal','maligno'),3)
AUC_svm_teL <- round(measureAUC(as.data.frame(pred_svmL)$prob.maligno, as.data.frame(pred_svmL)$truth, 'normal','maligno'),3)
AUC_svm_teP <- round(measureAUC(as.data.frame(pred_svmP)$prob.maligno, as.data.frame(pred_svmP)$truth, 'normal','maligno'),3)
AUC_svm_teR <- round(measureAUC(as.data.frame(pred_svmR)$prob.maligno, as.data.frame(pred_svmR)$truth, 'normal','maligno'),3)
AUC_values <- rbind(AUC_lg_te, AUC_svm_teL, AUC_svm_teP, AUC_svm_teR)
AUC_values <- as.data.frame(AUC_values)
AUC_values$Modelo <- c('Reg Log','SVM L','SVM P','SVM R')
colnames(AUC_values) <- c('Area debajo de la curva (AUC)','Modelo')
row.names(AUC_values) <- NULL
AUC_values <- AUC_values%>%dplyr::select(2,1)
# Imprimo resultados
kable(AUC_values)
| Modelo | Area debajo de la curva (AUC) |
|---|---|
| Reg Log | 0.889 |
| SVM L | 0.891 |
| SVM P | 0.895 |
| SVM R | 0.902 |
Con un AUC de 0.902 sobre el set de test, el Kernel radial es superior. Las diferencias en rendimientos, con esta configuración, no son demasiado grandes. El paso siguiente, sería comparar el rendimiento luego de una búsqueda exhaustiva de hiperparámetros para todos los modelos.
Específicamente, el “costo” principal de SVM tiene que ver con la pérdida de interpretabilidad directa.
Nota Bibliográfica Snippets de código fueron basados en: 1. https://www.r-bloggers.com/2019/10/support-vector-machines-with-the-mlr-package/ 2. https://www.datacamp.com/community/tutorials/support-vector-machines-r